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Abstract 

We study density waves in the flows of granular particles through vertical 

tubes and hoppers using both analytic methods and molecular dynamics (MD) 

simulations. We construct equations of motion for quasi one-dimensional 

systems. The equations, combined with the Bagnold's law for friction, are 

used to describe the time evolutions of the density and the velocity fields for 

narrow tubes and hoppers. The solutions of the equations can have two types 

of density waves, kinetic and dynamic. For tubes, we can show the existence 

of kinetic waves, and obtain the condition for dynamic waves for tubes from 

the equations. For hoppers, we obtain the solutions of the equations up to 

the first order of the opening angle, which also show the existence of kinetic 

waves. We reproduce density waves in the MD simulations for tubes. The 

waves are believed to be kinetic based on a few evidences, including a well 

defined flux-density curve. In MD simulations of flows in hoppers, we find 

density waves, which are also believed to be kinetic. 
05.40+j, 46.10+z, 62.20-x 
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I. INTRODUCTION 



Systems of granular particles (e.g. sands) exhibit many interesting phenomena 
The formations of a spontaneous heap 0-0 and convection cells |8|-12| under vibration, 
and the segregation of particles [13~n| are just a few examples. These phenomena are 
consequences of unusual, and often complex, dynamical responses of the system. Considering 
the complexity of the dynamics, one is tempted to first have comprehensive understanding 
of the dynamics of granular media in relatively simple geometries, and proceed to more 
complicated situations. Even for the simple geometries like shear cells, vertical tubes and 
hoppers, granular media still show complex dynamics. For example, granular particles in a 
shear cell show not only non-Newtonian behaviors [|19H 22[|, but also stick-slip motions [[23|] as 
well as density waves |24|.|25|]. Here, we study the flows of granular particles through vertical 
tubes and hoppers, whose geometries are as simple as shear cells, but are less well studied. 
In granular flows through vertical tubes, Poschel found that the particles do not flow 
uniformly, but form regions of high density which travel with velocity different from that of 
the center of mass. The mechanism for the traveling density pattern or the density wave 
is, however, not clearly understood. The density waves are also found in outflows through 
hoppers [[Z7],^|]. Especially, Baxter et al |28| show that the velocities of the density waves 
are dependent on the opening angle of hoppers, and even their directions can be changed. 
Furthermore, they find the density waves only when they mix some amount of rough sands 
with smooth sands. In MD simulations of hopper, the density field is found out to be 



non- uniform |29 |. Again, the mechanism for the density wave is not clearly understood. 
We first construct equations of motion for quasi one-dimensional systems. With the 



Bagnold's law for friction [19], the equations are used to describe the time evolutions of the 
density and the velocity fields for narrow tubes and hoppers. The solutions of the equations 
can have two types of density waves, kinetic and dynamic. The dynamic waves are similar 
to sound waves, while kinetic waves are due to the kinematics of the equations |p0[ . For 
tubes, we can show the existence of kinetic waves, and obtain the condition for dynamic 
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waves by a linear stability analysis. We also study the effects of static friction, and find 
that it changes the qualitative behaviors of the flows. For hoppers, we obtain the solutions 
of the equations up to the first order of the opening angle, and also show the existence of 
kinetic waves. Next, we study the systems using molecular dynamics (MD) simulations. For 
tubes, we reproduce density waves, and we find a few strong evidences that the waves are of 
kinetic nature, which includes a well defined flux-density curve and the density dependence 
of the velocities of the waves. We also study the existence of dynamic waves by changing 
the inelasticity of the particles, and do not find dynamics waves. In MD simulations of 
flows in hoppers, we find density waves, which are believed to be kinetic waves due to the 
strong correlation between the density and the velocity fields. For hoppers with periodic 
boundary conditions, we find additional density waves, which eventually dominates. We find 
a few evidences that the additional waves are of dynamic nature. We, however, do not find 
dynamic waves in hopper with open boundary conditions. The results obtained by these 
theoretical and numerical studies are compared with the experiments. 

This paper is organized as follows. In Sec. [D], we present theoretical results. Detailed 
discussions on the definitions and properties of dynamic and kinetic waves are given in Sec. 
|I1 A| . The kinetic and dynamic waves in tubes are discussed in Sec. |11B 1| and Sec. |11B 2 



respectively. Also, kinetic waves in hoppers are discussed in Sec. pi I C[ We present the 



numerical results in Sec. [IH| . First, we define the interactions between the particles in Sec. 
[ [I A| . The results for kinetic and dynamic waves in tubes are given in Sec. [I I IB], and the 



results for hoppers are discussed in Sec. |III C| . Brief summary as well as the limitations of 



this work are given in Sec. IV. 



II. THEORETICAL APPROACH 



3 



A. Dynamic and kinetic waves 



Here, we construct equations of motion for flows in one-dimension, which will later be 
used to describe the time evolutions of fields in vertical tubes and hoppers. Let us define 
p(x,t) and v(x,t) be the density and the velocity at position x and time t, respectively. If 
we require the mass of the system to be conserved, 

|p+|W=0. da) 

Other equation comes from the momentum conservation, 

d d 

P-QjV + P v g^ v = F{ x >t)> ( lb ) 

where F(x, t)dx is the total force acting on the mass within [x, x + dx] at time t. 

There are two distinct mechanisms by which density patterns can travel. One is due 
to the instability of uniform-density flows to density fluctuations ("dynamic wave"), and 
the other is of kinetic origin ( "kinetic wave" ) JJIJ . We now discuss both types of waves in 
detail, paying particular attention to their differences. In order to make the discussion more 
concrete, we use a specific form of the force 

F(x,t)=p-^-p-p^, (2) 

where /i and v are constants. The form of Eq. (|2]), motivated from a granular flow in a 
tube, is chosen such that the discussion becomes more clear. The followings are, however, 
applicable for any form of the force. 

We first consider the dynamic wave. The equations of motion Eqs. (|T|) with the force 
given by Eq. @ have a solution of uniform density p , 

p(x,t)= p , 

v(x,t)=p^. (3) 

We study the stability of the solution Eq. (|3]) by a linear stability analysis. For simplicity, 
we consider the case of p a = 1. The trial solution for the time evolution of a perturbation 
of Eq. (|) is 



p(x, t)= 1 + e p exp[i(kx — ut)] 
v(x, t) = 1 + e v exp[i(kx — ut)}, 



(4) 



where e p and e v are constants. We want to study the time evolution in a linear approximation, 
where we consider terms up to the first order of e. The analysis is valid only if a perturbation 
term is much smaller than unity. Substituting Eq. (|J) into Eqs. (|l|), and discard terms higher 
order than e, we obtain 



The density fluctuation travels with velocity Re(u)/k, where Re(u) is the real part of uj. 
If Im(uj) (the imaginary part of u) is positive, the perturbation grows with time, makes 
uniform- density flow unstable. This mechanism for creating density waves, which is based 
on the instability of the equations, is called dynamic waves. 

We now discuss the mechanism for kinetic waves. For an excellent introduction to the 
subject, please see Ref. P0| . Consider a case that the system is divided into two uniform- 
density flow regions with densities p a and p&, and the velocities of each region is determined 
by Eq. (|3|) (Fig. [I]). The equations of motion Eqs. ([I]) imply that density and velocity field 
(p(x,t) and v(x,t)) inside each domains do not change. The system evolve only by moving 
the interface. Let U be the velocity of the interface. We also set the initial position of the 
interface to be without losing generality. Then, the position of the interface at time t is 
Ut. At a given time, we choose an interval [— e, e] which includes the interface. We then 
integrate Eq. (pTaj) in the interval 



{uj — k)e p — ke v = 
[k + i(l - v)]e p + [(k — u) — ip]e v = 0, 



(5) 



whose solution is given by 



uj = 




(6) 
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(pa - p h )U + j b ~ja = 



(7) 
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where the flux j is defined to be pv. Therefore, the interface, or any density fluctuations, 
travels with velocity 

Pa - p b Ap 

Since the above mechanism for the travel of density fluctuations is a consequence of a kinetic 
equation (mass conservation, Eq. ([la])), it is called kinetic wave. Unlike dynamics waves, the 
fluctuations in kinetic waves can not be amplified, the system can only reorganize existing 
fluctuations. The fluctuations are created by initial conditions or by noises in the dynamics. 
Furthermore, kinetic waves, in most cases, decay into uniform flows |30 . 



In general, flows in 1-d have both kinetic and dynamic waves. If a system is stable under 
density fluctuations, existing dynamic waves will decay exponentially, and only kinetic waves 
can survive. On the other hand, if the system is unstable, the amplitudes of the dynamic 
waves will grow, and they eventually dominate over kinetic waves. In order to illustrate this 
point, we integrate numerically Eqs. ([[]) to get p(x, t) and v(x, t), where the force is given by 
Eq. (Q). The integration is done using the staggered leapfrog method [31] with the time step 



0.001. We also use a periodic boundary condition. The initial conditions are that p(x, 0) is 
a short pulse of higher density (1 %) on the top of background density of unity, and v (x, 0) 
is given by Eq. (§). We first study the case of p = 2, v = — 1. The behavior of the dynamic 
waves is studied by a linear stability analysis. Equation @ gives the dependence of uj on 
k to be 2k, —2i. There is one marginally stable mode of velocity 2, and one stable mode of 
velocity 0. On the other hand, the flux is given as 

j = pv = pfo^D/M (9) 

and the velocity of kinetic wave, determined by Eq. (Q), is approximately 2. Since the system 
is not unstable to density fluctuations, the kinetic waves, not greatly affected by the dynamic 
waves, will survive for a long time (Fig. 0(a)). We now study the case of p = —2, v — — 1. 
Following the same analysis, u(k) is determined to be 2k + 2i and 0. One unstable mode 
of velocity 2 and one marginally stable mode of velocity 0. The velocity of the kinetic wave 



is approximately 0. The dynamic wave of increasing amplitude (that of the velocity 2) is 
dominating over the kinetic wave (Fig. @(b)). 

B. Density waves in vertical tubes 

1. Kinetic waves in vertical tubes 

We consider vertical tubes of narrow width, which are filled with granular material. The 
particles will flow down due to gravity. This motion will induce the forces between the 
particles as well as the friction forces by the inner walls of the tubes. We assume that these 
forces follow the Bagnold's law The law, first found in a rapid granular flow, implies 
that the components of the stress tensor r„ are proportional to the square of the shear rate 

7, 

Tij = p B D 2 f i:j (p)^. (10) 

Here, pb is the density of the material which forms the particles, and p is the volume 
fraction of the particles. The density p is, therefore, ps times p. Also, D is the average 
diameter of the particles, fij(p) is a system dependent function. Bagnold predicted the 
shear rate dependence in rapid flows using a simple argument, which later confirmed by 
more elaborate calculations as well as experiments [0,[|^D]. In this paper, we will use the 
law except the case that the velocities of the particles are considered to be very small. 

Let x-axis (y-axis) be the direction parallel (perpendicular) to the tube, where the pos- 
itive x-direction is chosen to be upward. We consider the forces acting on the granular 
material contained in [x, x + dx}. Since we study only narrow tubes, we assume the material 
is homogeneous along the y direction. There is a gravitational force given by —pgWdx, 
where g is the gravitational acceleration, and W the width of the tube. There also is a 
friction by the wall, which is —r xy dx, and pressure by other particles —dr xx /dx Ddx. The 
total force F(x, t) becomes 
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F(x,t) = -p B pgW - S ign{v)p B D 2 f xy {p){^) 2 - D — [p B D 2 f xx {p){^) 2 ] 

d 

= -pBPgW - sign(v)p B f xy (p)v 2 - D—[p B f xx (p)v 2 ], (11) 

where sign(v) is the sign of v , and we assume the width of shear layer to be of order of D. 

Having found the force, we now proceed to study kinetic waves in the system. We first 
consider the case of no static friction. The steady-state velocity v s (p) of a uniform flow 
is obtained from the condition that the total force acting on the granular material is zero. 
Since we consider uniform density flows, the force free condition gives 



v 9 (p) = -yJpgW/fzyip). (12) 

The flux j(p) is given by p B pv s (p), then the velocity of a small density fluctuation U(p) 
becomes 

U(v) = ^ ~ ^ = iy^ ffxvW-pdfM/dP (13) 
Ap dp 2 s f xy (p) 



Therefore, the form of the force Eq. ( |lT|) allows density fluctuations to travel with velocity 
U(p). The condition for the existence of kinetic waves is the balance between the forces, 
and the waves are not very sensitively dependent on the exact form of the forces. Therefore, 
they are very likely to be seen in the experiments. 

We now consider the effects of static friction on the kinetic waves. By static friction, 
we mean one has to apply a finite force in order to break contacts between surfaces. There 
are two contributions to the normal stress on the wall. One is a static pressure P s , which is 
independent of the motion of the particles. The other is due to the collisions of the particles 
on the wall T yy , given by Eq. fllPl) . The total friction force is the minimum of pgW — r xy 
and p{j yy + P s ), where p is the friction coefficient. If the shear force (gravity and friction) 
is smaller than p times the normal force, the total force on the particles is zero, and the 



particles form a stagnant zone of no movement [0 . 

Depending on the parameters, there can be two behaviors of the system. If the gravita- 
tional force pgW is less than the friction due to static pressure pP s , the particles can not 
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move, and the system remains as a static column. If pgW > pP s , the particles will start 
move. The steady-state velocity v' s (p) for given packing fraction p is, 



v'sip) 



pgW-pP s /p B _ ^ 



\ fxyip) 

which is very similar to Eq. (|T2|). However, we have additional constraint that the shear 
force should exceed p times the normal force, which gives 

pgW - p B fxy(p)v 2 > pP s + pp B fyy(p)v 2 . (15) 

If we define v c (p) be the velocity that above equation becomes equality, 

V C (P) = V' s {p)— ) r . pr < V' s (p). (16) 

1 + Pfyyip)/f x y{p) 

Since the inequality Eq. dl5|) is not satisfied for v(p) > v c (p), the particles can not reach 
their steady state. Instead, they forms a stagnant zone. If we assume that the particles are 
stopped completely in the stagnant zone, the distances between the zones T>(p) are constant, 
which is the distance needed for the particles to be accelerated to v c (p) from the stopped 
state. We can calculate the distance T>(p) as follows. The velocity of the particles, obtained 
by integrating the Newton's equation with the friction force of Eq. flT5|), becomes 

uP 

v(t,p) = ^(p)tanhK(p)( 9 - ^)t\. (17) 

The distance covered by the particles during the interval [0,t c ] becomes T>(p), where t c is 
given by v(t c ,p) = v c (p). Integrating Eq. (|17D, we obtain 

^>- (g -^U in(i -^ (i8) 

where 7 = 1/(1 + pf yy (p)/ fxyip))- The distance T>(p) increases as p is decreased, and it 
diverges when p — > 0. 

We now summerize the results obtained in this section. If static friction is not present, 
the particles in vertical tubes form density waves travels with velocity U(p). When we take 
into account static friction, the system forms periodic stagnant zones with their separations 
V(p). 
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2. Dynamic waves in vertical tubes 



In this section, we study the stability of the uniform- density flows through vertical tubes 
under density and velocity fluctuations. We will follow the analysis presented in Sec. [11 A. 



We start with the equations of motion Eqs. (^) with the force given by Eq. (|TT|). Here, we 
do not consider static friction. The uniform- density flow solution with the density p Q is 

p(x,t)= p = p B p 



We study the stability of the solution by studying the time evolution of the perturbation 
from the uniform-density flow, given by 

p(x,t) = p + e p exp[i(kx — ut)} 

v(x, t)= v + e v exp[i(kx — ut)]. (20) 
Substituting Eqs. (P0|) into Eq. (|1|), and consider only up to the first terms of e, we obtain 

(kv - uj)t p + p Q ke v = 

Since the two equations in Eq. (ETJ) are valid for any values of e p and e v , the solution is 



_ i^L(f xy - tkDf xx )n + ik[g - ^(^k - tkD^)} = 0, (22) 
Po W KJ v J 1 [y W K dp dp n y ! 

where Q = uo — kv Q . The stability of the uniform solution is determined by the imaginary 
part of u. Unfortunately, since uj depends on the exact form of the unknow function fy, we 
can not determine the stability of the flow. Similar stability analysis on shear cells [^J and 



dissipative gases [33] shows that the system is unstable for small coefficient of restitution e. 



Especially, in ||33|| , this instability is traced back to the fact that pressure can decrease as 
the density is increased. Here, we can show that Eq. ( p2]) has an unstable mode if df xx /dp 
is sufficiently negative (smaller pressure for larger density), where the exact criterion is a 
complicated function of k, f xy and f xx . The stability of the flow is later checked by MD 
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simulations, where we will study the stability for various degrees of the inelasticity of the 
particles. 



C. Density waves in hoppers 

Consider a hopper with the opening angle 29, where the width at position x is given by 
W(x) = W + 2tan6 l x (Fig. 0(a)). Here, we consider only the positive ranges of x. The 
equations of motion for flows in hoppers are slightly different from those for tubes. Since 
the width W is dependent on x, the conservation of mass implies 

w §i p+ iL (Wfw) = o > (23) 



instead of Eq. (|la]) . Equation (lb) , which is a consequence of momentum conservation, still 



holds for hoppers, with a modified form of friction force. Since the sidewalls are tilted, the 
friction per unit length along the walls have two contributions. One is the component of r xy 
(friction force for tubes) parallel to the walls r xy cos 9. Forces on the walls by the internal 
pressure t xx give additional contribution t xx sin 9 cos 9, where sin 9 is due to the cross- section 
for the collisions of particles to the walls. (Fig. 0(b)). Therefore, the total force becomes, 

d 

F(x,t) = -pg + p B v 2 (f xy cos9 + f xx sin 9 cos 9) - p B D—(f xx v 2 ). (24) 

Since we assume that the system is homogeneous in the horizontal direction, we only consider 
hoppers of narrow width W . We also assume 9 <C 1, and the 9 dependences of physical 
quantities (density, velocity) are calculated up to the first order of 9. 

We now study the steady state properties of hoppers. Imposing the steady state condition 
on the mass conservation Eq. (|23|) gives 

JL(Wpv) = 0. (25a) 
ox 

Also, the momentum conservation Eq. ([Tb[) with the force Eq. (^4[), in a steady state, becomes 

9 1 D d 

VV ~dx V = ~~ P9 + w v2 (f x y cos ® + fxxSin9cos9) - — — (f xx v 2 ). (25b) 
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We want to know the density and velocity fields which satisfies Eq. fl25|). We start with a 
trial solution 



p(x)= p + A(x)8 

v (x)= v s + B(x)6, (26) 

where v s = — \Jp gW / f xy (p ), the steady state velocity for tubes Eq. (|l"2"D. Substituting Eq. 
(^6|) into Eq. fl25|) , and ignore terms of higher order than 8, 

TXr dA _„ dB 
v s W — + p W —= -2p v s 

KAJiXj \JbiAj 

Dv s^~^ + ( 2D VsU* + PoVsW )^£ + gW Q A - 2v s f xy B= v 2 J xx - 2xg Po . (27) 

Since f xx and f xy are also functions of p(x), one has to know the exact form of these 
functions to find the solutions of the equation Eq. (|27j). Here, we assume that f xx (p) can be 
approximated to f^ x p m for the range of densities found in a hopper. Similarly, f xy ~ f xy P- m 
Since 6 1, the changes of A(x) and B(x) to x will be also slow. We therefore consider 
variations up to the first order of x, that is, A(x) = A a x + A\ and B(x) = B Q x + B\. 
Substituting f xx ,f xy and A(x),B(x) into Eq. ([2TD, we obtain 

„PoV s 



v s A a + p B 



W 



2g Po = {mvlr xyV T X ~ gW )A + 2v w f° yP ™B 
v sPo W B = -gW A l + v 2 s p^f° xx + mf xy v 2 s p n - 1 A 1 

+2f° xy v sP ™B 1 - Df° xx v s p™{2B + ^^). (28) 

" Po 



There are only three conditions in Eq. (|28|) for four unknown variables. The other condition 
comes from the freedom in choosing the origin of the x axis. We set the origin so that 
B 1 = 0. Then, the solution of Eq. @ is 

„ 2m v x 



3-m W 



W 3-m 

a '=(^4^- [ 3^ (p °- / ^ ) - /> ™ 1 ' (29) 
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Intuitive solutions for hoppers are that the magnitude of the velocity increases, and the 
density decreases as x is increased. Although the solution Eq. ( |29| ) becomes intuitive one 
for < m < 3, it can also have very different behaviors for other ranges of m. 

Having obtained the steady state for hoppers, we now proceed to study kinetic waves. 
The density and velocity fields are described by Eq. with one free parameter p Q . Consider 
two regions of density, where p Q is chosen to be p\ and P2, respectively. The velocity of the 
interface Uh between the regions, which can be calculated following the way described in 



Sec. II A 



is 



jj = Aj = A\piV s (pi) + VsPiAjjpije + (piB (pi) + v s (p, i )A o (p i ))x0} 
h ~ Ap " &\piV s {pi) + (A ( Pi ) + A l {p i )x)6} * 1 ' 

We consider small fluctuations of density: p\= p and Sp = p 2 — p <C 1. The velocity, up to 
the first order of 8, is 

+^-( P B (p) + v s {p)A o (p))x0 - —(pv s {p))—{A {p)x + A 1 (p))9. (31) 
dp dp dp 

The above velocity Uh is calculated at position x. One should note that the density at the 

position is not p but p + (A (p)x + Ai(p))0. We now compare the velocity of density waves 

in hoppers and tubes. The velocity of density waves in tubes U t at the above density is 

d f dv 8 {p) dv s (p) d 2 v s (p). 

which is give by Eq. (|13|). 

Therefore, density waves also exist in hoppers, and their velocity is given by Uh{0,p, x) 
Ut + C(p, x)9, where C(p, x) is a complicated function of x and p. 



u t = i:(pvs(p)) + + p^^ +P^f l )(Mp> + Mp)% (32) 



III. MOLECULAR DYNAMICS SIMULATIONS 

A. Interactions between particles 

We discuss the interaction between the particles used in the MD simulations of granular 
flows. The force between two particles i and j, in contact with each other, is the following. 
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— * — * — * — * 

Let the coordinate of the center of particle i (j) to be Ri (Rj), and r = Ri — Rj. In two 
dimensions, we use a new coordinate system defined by the two vectors n (normal) and 
s (shear). Here, h = r/\r\, and s is defined as rotating n clockwise by tt/2. The normal 
component F^_^ of the force acting on particle i by j is 

= k n (cii + aj - \r\) 3/2 - j n m e (v ■ n), (33a) 

where Oj (a,) is the radius of particle i (j), (rrij) the mass of particle i (j), and v = df/dt. 
The first term is the Hertzian elastic force, where k n is the elastic constant of the material. 
And, the constant 7„ of the second term is the friction coefficient of a velocity dependent 
damping term, m e is the effective mass, rriimj/(mi + rrij). The shear component F?_^ is 
given by 

= -7 s m e ({7- s) - sign(5s) mm(k s \Ss\, nlF^l), (33b) 

where the first term is a velocity dependent damping term similar to that of Eq. ( |33a| ). The 
second term is to simulate static friction, which requires a finite amount of force (/i-Fjl^) 



to break a contact | 34| . Here, \x is the friction coefficient, 5s the total shear displacement 
during a contact, and k s the elastic constant of a virtual spring. There are several studies 
on granular systems using similar interactions. However, only a few of them [j34| -|3"6| include 
static friction. A particle can also interact with a wall. The force on particle i, in contact 
with a wall, is given by Eqs. ([33|) with cij = oo and m e = mi. A wall is assumed to be 
rigid, i.e. it is not affected by the collisions with particles. Also, the system is under a 
gravitational field g. We do not include the rotation of the particles in present simulation. 



A detailed explanation of the interaction is given elsewhere |36[ 



B. Density waves in vertical tubes 

1. Waves without static friction 



We first simulate granular flows without static friction. Thus, we set /x = 0, and the 
shear force is only due to the velocity dependent friction term in Eq. ( |33b|) . We study the 
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system in two dimensions. Tubes are modeled by two parallel sidewalls of length L, and 
separated by distance W, between which particles flow (Fig. We use a periodic boundary 
condition in the vertical direction. Particles come out of the bottom of the tube are fed into 
the top. In order to avoid a haxagonal packing formed by monodisperse particles, we use 
polydisperse particles, whose radii are drawn from the gaussian distribution of mean 0.1 and 
width 0.02. We initially arrange particles to be equally spaced along the vertical direction, 
and calculate the positions and the velocities of the particles at subsequent steps using a 
fifth-order predictor-corrector method. 

In Fig. [5], we show the time evolutions of the densities and the velocities of the particles 
with W = 1, L = 15 and the number of particles N = 225. The density plots are made 
as follows. We divide the tube into several regions (bins) of equal height (typically, 5 times 
particle diameter), and count the number of particles in bin %. We set the grayscale of 
each bin to be proportional to rij. We choose white for rii = di, and black for rij = d u , 
where d\ (d u ) is the lower (upper) bound for nj. If rii is smaller than di or larger than d u , 
the grayscale is chosen to be white or black, respectively. In Fig. |5|(a), we use d\ = and 
d u = 30. The density field at a given time step is plotted as one horizontal line, where boxes 
of different grayscale represents bins of the tube. Here, the leftmost box corresponds to the 
bin at the bottom of the tube. The velocity plot is made using the same procedure as above, 
except the grayscale is proportional to the average vertical velocity V{ in bin i. In order to 
enhance the contrast, we subtract the center of mass velocity from t>j. We choose white for 
Vi and black for —v\. Here, we set v\ = 60. The time step is chosen to be 5.0 x 10. ~ 5 The 
time interval between the successive raws in the density and velocity plots is 100 iterations. 
The parameters for the simulation are k n = 1.0 x 10 6 , k s = 1.0 x 10 4 and 7„ = 7 S = 500 
between the particles. Between the particles and the walls, we use k n = 5.0 x 10 6 in order 
not to allow the particles go through the walls, while the other parameters are kept to be 
the same. 

In Fig. [|, one can see a region of high density is being formed from the homogeneous 
system |$7|. Also, a high density region may be split into two, or two regions may merge to 
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form one region. However, for the most of time, these density fluctuations just travel with 
almost constant velocity. These traveling density patterns are first observed in simulations 
by Poschel pB| . Comparing Fig. |5](a) and (b), one can notice correlations between the 
density and the velocity V{. The particles seem to travel slower (faster) in high (low) 
density regions, which is very similar to traffic flows. This correlation is one hint that the 
density waves may be of kinetic nature. In order to systematically study the correlation, 
we measure the density dependence of the flux in a steady state. We choose a bin, and we 
calculate the friction and gravity force acting on the bin. Since we want to measure the flux 
in a steady state (the total force is zero), we require the total force (sum of the friction and 
the gravity) is smaller than r times the gravity, and we discard bins if the requirement is 
not fulfilled. For bin i in a steady state, we measure the total flux, defined as ji = ViUi. We 
calculate the density dependence of the average flux ji for a system of fixed total number 
of particles N. Here, the averages are taken over time. We measure the flux-density curve 
for several different values of r ranging from 0.1 to oo. The results are not very sensitive 
to the values of r, when we study the systems in a steady state. We set r = 1.0 from now 
on. We also measure the curve for several values of N = 150, 225, 280, 337. For large N, we 
have accurate estimate of ji for large values of rij, but poor one for small ?v The situation 
is opposite for the systems of small N. For the intermediate values of rij, however, all the 
systems give good estimates, which agree with each other. In Fig. |6|, we show the flux-density 
curve averaged over the four values of N. The fact that we have a well defined flux-density 
curve suggest the system is in a steady state, which implies that kinetic wave is sufficient 
for the description of the evolutions of the system. Furthermore, the curve resembles that 



of a traffic flow, which is considered to be one of the typical examples of kinetic waves |30 

One additional evidence that the density waves are of kinetic nature is the dependence of 
the velocity of the waves on the average density. The velocities of kinetic waves are given by 
Eq. (|§p, which are the slope of the flux-density curve, for small density fluctuations. From 
Fig. |6], we expect that the velocity is a large negative value for a small density, approaches to 
zero, and becomes a large positive number as the density is increased. We directly measure 
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the velocities of the waves from the slopes of high density regions like the one shown in Fig. |5|. 
In Table. |, we show the average velocities for several values of N, where the average density 
< rii > is given by N/L. In the table, one note the velocity is negative (—41) for small N, 
and is increased to positive (113) for large N, which is exactly the way predicted by the 
theory of kinetic wave. Furthermore, the measured velocities are consistent with the local 
slopes of the flux-density curve Fig. || although the slopes can not be accurately determined 
due to the large error bars. Based on the above evidences, as well as the theoretical argument 



given in Sec. [II B 1| , we conclude the above density waves found in vertical tubes are kinetic 
waves. 

We also want to discuss the origin of fluctuations in the system. As shown above, 
density fluctuations (waves) are formed from an uniform density system. But kinetic waves, 



as discussed in Sec. |I A| , can not create fluctuations. Also, some of the waves are split into 



two waves, which can not be described by the evolutions of kinetic waves alone. There must 
be some sources of fluctuations or "noises" in the system. Since the system is deterministic, 
one might think the system can not have noises. The "noises" come from the fact that the 
equations of motion Eq. ([!]) as well as the form of the friction force are relations between 
averaged quantities. The fluctuations around their averaged values, especially relevant in 
small scale descriptions of systems, are identified as "noises" . 

2. Waves with static friction 

We next study the flows through vertical tubes with static friction. As discussed in 
Sec. [II B 1| , we expect two types of behaviors depending on the friction coefficient /i. If 
fi > pgW/P s , the system can not move, and stays as a static column. Otherwise, the particles 
start move, and forms periodic stagnant zones whose separation T>(p) is determined by Eq. 



18. The setup we used in the MD simulations is exactly the same as the previous one (Fig. 
^|), where static friction is introduced by choosing a non-zero fi. Between the particles, we 
use upp = 0.5. We also set the shear friction coefficient 7 S to be zero. We do simulations 
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with W = 1 and L = 15, starting with iV = 225 particles. All the other parameters are kept 
to be the same as above. In the simulations, we find three types of behaviors depending 
on the friction coefficient between the walls and the particles fiwp'- (1) If I^wp > 1-0, the 
particles can move initially, but they eventually form static column(s). (2) On the other 
hand, if \i-wp < 0.6, the particles constantly increase their speeds, and do not go into a 
steady state until the end of simulations (50 000 iterations). (3) In the intermediate regime, 
0.6 < Hwp < 1-0, we find steady states, where density waves travels with almost constant 
velocities. The static structures found in the first regime are precisely what is expected 
from the theory. In the second regime, however, the systems do not reach steady states, in 
contrast to the theory. One possibility is that the time needed to reach a steady state is 
larger than the simulation time. To check this possibility we do longer simulations of the 
system with fiwp — 0.6, which does not reach a steady state in 50000 iterations. We find 
the system does reach a steady state at around 100000 iterations. 

In the third regime, the simulations with L = 15 do not show the expected periodic 
stagnant zones but traveling density fluctuations. For example, we show the density evolu- 
tions in a tube of L = 15 in Fig. 0(a), where fiwp — 0.9 and the time interval between the 
successive raws are 0.0025. A region of high density travels with almost constant velocity of 
—86 ± 13. There are few possibilities to understand the discrepancy. First of all, there can 
be a problem of commensurability. In general, the length of the tube L is not an integer 
multiple of the stagnant zone separation V{p). Thus, the distance between the clogged zone 
can not all be T>(p) due to the periodic boundary condition. It can be shown, following the 
argument in Sec. [II B 1| , that the distances between the stagnant zones are all T>{p) except 



one, which is smaller than T>(p). This configuration, however, can not be a steady state. 
The stagnant zone just below the small separation becomes unstable, since the incoming 
velocity of the particles are smaller than v c (p). These particles travel further down to reach 
v c (p), then form a stagnant zone. In effect, the stagnant zones travels down, if L is not 
an integer multiple of T>(p). There also can be problems of "noises." As discussed before, 
there are noises to the continuum description we have been using. Due to these noises, the 



separations between the zone are not all equal, but fluctuate around T>(p). Furthermore, 
the noise also cause the stagnant zones to move around. Considering the above possibilities, 
we expect that several stagnant zones in a long tube drift downwards, whose separations 
fluctuate around T>(p). In order to check the possibilities, we study density waves for several 
values of L, keeping the average density constant. The evolution of density field for L = 30 
is shown in Fig. 0(b). One can see a region of high density is traveling downwards with 
the velocity around —25. The results for both L = 15 and L = 30 are consistent with the 
prediction. Since there are only one region of high density, we expect T>(p) is larger than L. 
We then expect the velocities of the density waves are always negative, and their magnitude 
decreases as one increases L, thereby approaching T>{p). The situation for L = 45 is quite 
different, as shown in Fig. |7|(c). The high density region seems to travel upwards, and the 
velocity seems to fluctuate more. These two features are also found in a few additional runs 
we study for L = 45. The fluctuations can be caused by the "noise," but the trend of moving 
upwards can not be explained. We also studied the system with fj, WP = 0.7 and 0.8, and 
find essentially the same. 



3. Comparison with the experiments 

We now discuss the results of the experiments with vertical tubes, and compare them with 
the theoretical and simulational results obtained above. We first discuss the experiments. 



Density waves in vertical tubes are first found in the experiments by Poschel |26] . He found 
two types of waves: (1) Regions of large densities occur at random positions, and they 
travel with non-zero velocity. (2) The separations between the regions of large densities 



seem to be about the same, and they fluctuate around certain positions |38]] . The conditions 
needed to obtain each types of waves are unknown. The flows are also studied in vacuum, 
and surprisingly, the density waves disappear p9| , ^0| . It is still not understood how air 



affects the formation of density waves. Here, we discuss one of the possible mechanisms. 
The conditions for forming kinetic waves, as discussed in Sec. [11 B 1| , are fairly simple. One 
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important condition is that the friction force should balance the driving force, which in this 
case, gravity. It is possible that the friction force by the sidewalls is too small to balance 
gravity. Consider a block of particles falling down in a tube. If the density of the block is 
large enough, air can not easily pass through the block, and the air pressure just behind 
the block can be smaller than that of the front. The pressure difference gives rise a force to 
slow down the block, which acts as an additional friction force. Since the pressure difference 
is expect to increase by increasing the density and increasing the velocity of the block, the 
friction force by air can balance gravity at high velocity, and the balance produces kinetic 
waves. As one can see, this argument is largely speculative, and should be checked by careful 
experiments. Especially, it should be check whether there is pressure difference between the 
front and the behind of a moving granular blocks of high density, and whether the difference 
is enough to balance the gravity. 

Assuming the above mechanism to be hold, how can one understand the two types of 
waves found in the experiments? The two waves are readily compared with the predicted 
kinetic waves with and without static friction. If V{p) is larger than the length of the tube, 
we do not expect to see any stagnant zone (for open boundary condition). Therefore, we 
effectively see the system without static friction, which produces traveling kinetic waves. 
However, if T>(p) is much smaller than the tube length, we do see periodic stagnant zones. 
The distance V(p) is a decreasing function of the friction (both by the sidewalls and air), 
and it is possible that the tube which produces the periodic zones have larger friction than 
others. This possibility again is speculative, and should also be checked by more controlled 
experiments. For example, the friction can be changed by using the sidewalls of varying 
surface roughness. If air is responsible for an additional friction, one can also change the 
friction by controlling the pressure of air. 
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4- Dynamic waves 



We conclude this section by discussing dynamic waves in granular flows through tubes. 
As discussed in Sec. LIB 2j , we expect the granular flows become less stable under density 



fluctuations for smaller coefficient of restitution e. Since e is small for large normal damping 
coefficient j n , we expect the system to be less stable for large 7 n . We simulate the system 
with several values of 7„ = 5 x 10 2 , 1 x 10 3 , 2 x 10 3 , 3 x 10 3 , and study the stability. In 
Fig. ^, we show the density fields of tube with L — 15 and N = 225 for several values of 
7 n . Here, the time intervals between the successive raws are 0.0025. We do not find any 
sign of dynamic waves, even for longer simulations. We also repeat the simulations with 
higher density iV = 337, and find that the results are essentially the same. The coefficient of 
restitution e depends on the relative velocity between the colliding particles for the Hertzian 
contact force. Considering the scale of velocities found in tubes, e is roughly estimated to 
be 0.1 for j n = 3 x 10. 3 It is possible that the dynamic instability occur at even larger values 
of 7 n . The limitation of the constant timestep type algorithm we are using is that we have 
to decrease the timestep to prevent a numerical instability, if we want to increase 7„. For 
example, the timestep is chosen to be 5 x 10~ 5 for 7 n = 5 x 10 2 , and 1 x 10 -5 for 7 n = 3 x 10 3 . 
Therefore, we can not simulate systems of arbitrarily large j n . One have to use event driven 
type algorithms (for example, Ref. [[HJ]) to overcome this limitation. 

C. Density waves in hoppers 



In Sec. [II C| , we presented a theory which predicts the existence of kinetic waves in 
granular flows through hoppers. In that section, the velocity of the kinetic wave in hoppers 
of opening angle 29 is shown, up to the first order of 9, to be Uh(9,p, x) = U t {p) + C(p, x)9. 
Here, p is the packing fraction and x the position, and U t (p) the velocity of kinetic waves in 
tubes, and C(p, x) is a complicated function of p and x. We thus expect Uh for a hopper of 
small 9 is not very different from that of tubes. The conditions for the existence of dynamic 



21 



waves are too difficult to be obtained from the theory. 

We study hoppers of length L = 15, bottom width W = 1.0, and several values of the 
opening angle 28. We apply a periodic boundary condition in the vertical direction. The 
particles come out of the bottom are again fed into the top. The boundary condition, which 
is not natural for the hopper geometry, can introduce some artifacts to the system. We later 
check the results by comparing with those from the open boundary condition. The main 
reason for using the periodic boundary condition is that we can simulate the systems for 
longer time, which in effect study the systems of larger sizes without actually increasing the 
length L (and the number of particles). We initially arrange particles as a square lattice in a 
hopper. The lattice constant is 0.1, the average radius of the particles. The average number 
of particles per unit area is very close to 25, the maximum density allowed for the square 
packing. The randomness in the initial configuration is introduced by the polydispersity of 
the particles. The interaction parameters are chosen to be exactly the same as those for 
tubes. 

We first study the system without static friction. The static friction coefficient /i is 
chosen to be 0, and the velocity dependent shear friction term ^ s to be 500. In Fig. |^, we 
show the time evolutions of the density and the velocity fields in hoppers. The plots are 
made in the same way as Fig. ||] except the fact that the density is obtained by dividing the 
number of particles in a bin by the area of the bin, which is not constant for a hopper. The 
grayscales are chosen to show clear contrast between the low and high density (velocity) 
regions by controlling d\ and d u (v[). 

In Fig. |9|(a) and (b), we show the density and the velocity fields for a tube (8 = 0) in 
order to serve as a reference to compare with those for hoppers. The time interval between 
the successive lines are 0.005. In Fig. |9](c) and (d), the fields for a small angle hopper 
(8 = 1°) are shown. In the density plot (Fig. |9|(c)), there are two waves in the hopper, one 
travels upwards and the other downwards. Comparing with the velocity plot (Fig. |](d)), the 
densities of the upward wave are strongly correlated to the velocities, while there seems to 
be no such correlation for the downward waves. The correlations found in the upward waves 
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are qualitatively the same as that of kinetic wave, i.e., the particles travel slow (fast) in the 
regions of high (low) density. It is quite possible that the upward wave is of kinetic nature. 
Also, the velocities of kinetic waves in hoppers with small opening angle are expected to be 
very close to that of tubes. The fact that the velocity of the upward waves (113 ±4) are very 
close to that of the tube (93 ± 10) is one other support that the wave is kinetic JIT]. We now 
consider the waves travel downwards. The waves, not only show no correlations between the 
density and the velocity fields, but also travels in the opposite direction to the kinetic waves 
in the tube (Fig. H(a)), which suggests that the waves are not probably kinetic. Furthermore, 
the density contrasts of the waves are constantly increasing, which is not possible for kinetic 
waves. The increments can be more clearly seen in hoppers of larger opening angle as 
shown in (e)-(f) (9 = 4°). The downward waves initially coexist with the kinetic waves, 
but eventually dominate the system. The above facts are still true for the largest 9 = 10° 
we study. The properties of the downward waves listed above suggest that the downward 
waves are of dynamic nature. According to the argument used for the dynamic waves in 
tubes, dynamic waves will be more easily formed using non-elastic particles. We check this 
possibility by studying a hopper of 9 = 1° for different values of j n . In Fig. [LT], we show the 
density fields for 7„ = 1 x 10 3 and 2 x 10 3 . Comparing with that of 7 = 5 x 10 2 (Fig. |S|)(c), 
one can see the waves for larger 7^ are indeed formed earlier, and their intensities (density 
contrast) are larger. This add one more support that the downward waves are dynamic. We 
also measure the velocities of the dynamic waves for several values of 9. As shown in Table 
HI the magnitude of the velocity is decreased as 9 is increased. Following the trend, it is 
quite possible that the velocity becomes zero for finite 9, and change its sign. 

We now study hoppers with the open boundary condition. In Fig. |ll|(a) and (b), we show 
the density and the velocity fields for a hopper of 9 = 1° and L = 30. Here, we turn off static 
friction (/i = 0), and we set 7 S = 5 x 10 2 . One can hardly see any fluctuations in the density 
field, but one can see some traveling patterns in the velocity field. We now turn on static 
friction. We set 7^ = and fi 7^ 0. The density and the velocity fields with /i = 0.5 are shown 
in Fig. |ll|(c) and (d). One can now see upward traveling density waves. Furthermore, the 
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density field have strong correlations with the velocity field, suggesting that the waves are 
kinetic. Even though kinetic waves are present in the system independent of static friction, 
the amplitudes of the waves are too small to be visible without it. Static friction provides 
an effective mechanism for creating large density fluctuations. Also, the density waves are 
not present in the hoppers of small L. The absence of density waves can be caused by the 
fact that the particles need to travel certain distances to reach a steady state around which 
the density waves can only be seen. All these results remain valid for larger values of 9 we 
have studied. For example, we show the case of 9 = 5° in Fig. |TTKe)-(f). 

We also search for the dynamic waves in hoppers with the open boundary condition. We 
simulate the system with several values of j n , and check the instability of creating dynamic 
waves. For the range of 7„ we have studied (5 x 10 2 ~ 3 x 10 3 ), we do not find the instability. 
Again, the simulations for larger values of 7 S are limited due to the algorithm we are using. 
Even though we do not rule out dynamic waves in hoppers in general, the dynamic instability 
found above seems to be an artifact of the periodic boundary condition. 

We now compare the results with the experiments. Baxter et al found density waves in 
a hopper only if a certain amount of rough sands is mixed with smooth sands [J28| . The 
role of the rough sands is not clearly established. It is possible that the rough sands forms 
local "arches" , thereby increasing density fluctuations. In that case, the role played by the 
rough sands is the same as that of static friction, which provides density fluctuations to 
maintain kinetic waves. Also this is consistent with the fact that the amplitudes of the 
density waves are a smooth function of the fraction of the rough sands, which suggests 
that the density waves in the experiments are kinetic. The suggestion, of course, has to 
be checked by experiments, for example, by measuring the correlations between the density 
and the velocity fields. 
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IV. DISCUSSION 



We have presented theoretical and numerical evidences that the density waves found in 
the simulations of Poschel is of kinetic nature. However, the density waves found in the 
experiments are not fully understood. The first and foremost problem is to find the form of 
the friction force. In the MD simulations, the friction force is generated by the collisions of 
particles to the sidewalls. In experiments, the collisional friction force seems to be too small 
compared to gravity, and another friction force related to air gives dominant contribution. 
We proposed a mechanism how air can generates a friction force, which should be checked by 
experiments. We want to emphasize that since the existence of kinetic wave is not strongly 
dependent on the details of the friction force, the density waves found in the experiments 
are very likely to be kinetic waves. 

We find conditions for the existence of dynamic waves in tube. The conditions depend 
on details of the force, which we do not know. We numerically search for dynamic waves 
by increasing j n , since we expect dynamic waves can be more easily created for inelastic 
particles. Even for the largest 7„ we studied, we can not find the dynamic waves. It is 
possible that dynamic waves occur for even larger values of j n . We are not able to check 
the possibility, since the algorithm we are currently using becomes very inefficient for larger 
values of j n . 

We also have present theoretical and numerical evidences that there are kinetic waves 
in hoppers. Especially, the kinetic waves in hoppers with the open boundary condition are 
visible only with static friction. This can be readily compared with the fact that one need 
finite fraction of rough sands to observe the density waves in the experiments. Here, the 
role of rough sands can be creating large density fluctuations to maintain kinetic waves just 
like static friction. The suggestion that the density waves in the experiments of hoppers 
are kinetic waves should be checked, for example, by studying the correlation between the 
density and the velocity fields. 

We are not able to do a linear stability analysis for hoppers, due to the complicated 
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density and velocity fields in the steady state. In the simulations with the peridic boundary 
condition, we find another density waves, which we believe to be dynamic on a few evidences. 
Since we do not find the dynamic waves with the open boundary conditions, we think the 
above dynamic wave is an artifact of the boundary condition. However, we do not rule out 
the dynamic waves in hoppers, especially for large values of j n . 

I thank Michael Leibig, Hans Herrmann and Thorsten Poschel for many useful discus- 
sions. 
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FIGURES 



FIG. 1. The time evolution of the interface between the uniform-flows with density p a and pb- 

FIG. 2. Time evolution of density fluctuation with the force given by Eq. (2): (a) Dynamic 
waves are decaying, while the kinetic waves survives (// = 2,v = —1), (b) Dynamic waves are 
growing, and take over kinetic waves (// = — 2,v = — 1). 

FIG. 3. (a) Hopper with the opening angle 29. (b) The friction force per unit length on the 
wall is sum of T xy cos9 (left) and t xx cos 9 sin 9 (right). 

FIG. 4. Vertical tubes in two dimension with width W and length L under gravity. We apply 
a periodic boundary condition in the vertical direction. 

FIG. 5. Time evolutions of (a) density and (b) velocity fields of particles in a tube of 
W = 1,H = 15. Fields at a given time is shown as a horizontal line of boxes. The grayscale 
of each box is proportional to the density and the velocity in that region of the tube. Regions of 
high density are formed, and travel with almost constant velocity. 

FIG. 6. The flux-density curve for a tube of W = 1, H = 15 averaged over time and different 
values of N. The parabolic shaped curve resembles that found in a traffic jam. 

FIG. 7. The evolution of density fields in tubes with p = 0.5 and (a) L = 15, (b) 30 and (c) 45. 
High density regions travel downwards in (a) and (b), but travel upwards with larger fluctuations 
in (c). 

FIG. 8. Density fields for tube with L = 15 and N = 225 for several values of 7 n : (a) 
7„ = 1 x 10 3 , (b) 2 x 10 3 and (c) 3 x 10 3 . There is no sign of dynamic instability. 

FIG. 9. The density and velocity fields for hoppers with different opening angle, (a)-(b) 9 = 0°, 
(c)-(d) 1°, (e)-(f) 4°. For small 9, there are two waves traveling in opposite directions, where the 
downward waves eventually dominate. 
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FIG. 10. The density fields for a hopper of 9 = 1° with (a) 7n = 1 x 10 3 and (b) 2 x 10 3 . The 
intensities (density contrast) of waves are larger for larger j n . 

FIG. 11. The density and velocity fields for a hopper with the open boundary condition. In 
(a)-(b) the fields for a hopper of 8 = 1° without static frcition, (c)-(d) for 9 = 1° with static 
friction, and (e)-(f) for 9 = 5° with static friction are shown. 
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TABLES 



TABLE I. The velocity of kinetic waves in a tube for several values of the average density 



N <rii> Velocity of kinetic wave 

150 10.0 -41.0 ± 2.0 

225 15.0 5.0 ± 9.0 

280 18.7 12.0 ± 11.0 

337 22.5 113.0 ± 4.0 



TABLE II. The velocity of dynamic waves in a hopper for several values of 8. 

8 (in degrees) Velocity of dynamic wave 

1.0 -49 ± 3 

2.0 -40 ± 5 

4.0 -26 ± 4 

6.0 -21 ± 5 

8.0 -18 ± 4 
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